---
title: "Part 4: Deconstruction Wood Sensitivity Analysis"
subtitle: "2018 Deconstruction Data Net Impact Analysis using OSU average softwood density factor"
author:
- Andey Nunes, MS
- Research Analyst 2
- Oregon DEQ
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  html_document:
    df_print: paged
    toc: yes
    toc_depth: '3'
  word_document:
    toc: yes
    toc_depth: '3'
---

# Introduction, Project Goals, and Notebook Setup

The 2018 Deconstruction Data Analysis uses prepared and cleaned project data from residential single family homes removed under a City of Portland Deconstruction permit. The goal of the data analysis is to quantify the net environmental benefits resulting from avoided disposal of materials due to salvage/reuse (measured as Global Warming Potential and Primary Energy Demand impacts). This document is part four of the project where previously constructed data sets are used to determine and graphically summarize impacts under the deconstruction and demolition scenario and the net benefits of deconstructing residential homes over the hypothetical demolition of those structures.

The analysis is intended to be reproducible using R coding procedures and employs the following R Packages:
```{r packages, echo=F}
# create a list of packages to be installed and active (required) for use in the notebook and resulting reports
packages <- c("corrplot", "ggpubr", "ggrepel", "ggthemes", "grDevices", "knitr", "tidyverse")
# print this list as a table
knitr::kable(packages, col.names = "package name")
```

Prior to executing code and producing a report, the packages must be installed and accessed, and basic features of the document set up by specifying code chunk defaults.
```{r document setup, warning = F, echo=FALSE}
# set defalut code chunks to "echo=TRUE" to display the code chunk inline with the narrative text; and report numbers to a default of three significant digits.
knitr::opts_chunk$set(echo = TRUE, fig.width = 8, fig.height = 5)
options(digits = 1)
options(scipen = 999)

lapply(packages, require, character.only = T)
options(xtable.comment = F)

# DEQ color palette for graphics using approximate match colors from grDevices package

DEQ_pal <- c('aquamarine4', 'steelblue4', 'lightseagreen', 'yellowgreen', 'darkorange1',
            'darkseagreen3', 'slateblue1', 'powderblue', 'khaki', 'lightsalmon',
            'seagreen4', 'deepskyblue4', 'darkslateblue', 'magenta4', 'palegreen4',
            'cyan4', 'goldenrod2', 'indianred3', 'seagreen')


```

Specify parameters for saving graphics.
```{r graphics parameters}
# designate a folder name for saving graphics and the default size specification
graph_path <- file.path("graphs/")
width <- 8
height <- 5
```


## Data Import & consolidation 

Previously developed data sets are imported and cleaned/manipulated as needed, the index values that R assigns when saving *csv* files is dropped, units of the values are formatted to be consistent across all data tables, and contractor names are replaced with a random number for anonymity in the final report.

For worker and equipment impacts that are associated with projects but not a specific material category (which is what all observations in our `decon` and `demo` data sets are up to this point), we have to assign compatible names within the existing variable categories in order to create observations of each per-project based transport impact (i.e. workers and equipment). For this the `worker_equipment_impacts` data values will be assigned the `site_transport_impact` variable name, since job-site related transportion based impacts is the best fit description for the inclusion of these values. Additional names are given in other variables in order to distinguish these values from the material impacts and material related transport impacts (including site and EOL based transport). The `worker_equipment_impacts` table will have the following variables and assigned variable states (values) added prior to being combined with the `decon` and `demo` scenario data frames:

| variable name         | value (assigned variable state)   |
|-----------------------|-----------------------------------|
| `decon_equipment`     |  0                                |
| `SimpleEOLname`       | *per_project_transport*           |
| `disposition`         | *none*                            |
| `quantity`            |  1                                |
| `quantity_units`      | *per_project*                     |
| `material_impacts`    |  0                                |
| `EOL_transport_impact`|  0                                |


```{r data cleaning, warning = F}
# cannot use read_csv for reading in these files because of parsing error on some of the EOL transport dropbox values
decon <- read.csv("intermediary/final_deconstruction_scenario_with_transport.csv", stringsAsFactors = F)
demo <- read.csv("intermediary/final_demolition_scenario_with_transport.csv", stringsAsFactors = F)
worker_equipment_impacts <- read.csv("intermediary/per_project_worker_equipment_transport_impacts.csv", stringsAsFactors = F)
house_weight <- read.csv("intermediary/house_weight.csv", stringsAsFactors = F)
house_weight$project <- as.character(house_weight$project)
house_weight$contractor <- as.character(house_weight$contractor)

# standardize impact_units
decon$impact_units <- str_replace(decon$impact_units, "kgCO2e", "kg CO2e")

# clear index columns and add the activity scenario name
decon <- decon[,-1] %>%
   add_column(scenario = "decon", .before = 1)

demo <- demo[,-1] %>%
   add_column(scenario = "demo", .before = 1)

# adjust worker_equipment_impacts to fit the structure of the decon and demo activity scenario data frames
worker_equipment_impacts <- worker_equipment_impacts %>%
   add_column(decon_equipment = 0, .before = 7) %>%
   add_column(SimpleEOLname = "per_project_transport") %>%
   add_column(disposition = "none") %>%
   add_column(quantity = 1) %>%
   add_column(quantity_units = "per_project") %>%
   add_column(material_impacts = 0) %>%
   add_column(EOL_transport_impact = 0) %>%
   gather(`decon_workers`, `decon_equipment`, `demo_workers`, `demo_equipment`, key = "scenario_deconMaterialName", value = "site_transport_impact") %>%
   separate(scenario_deconMaterialName, into = c("scenario", "deconMaterialName"), sep = "_") %>%
   select(c("scenario", "house_age", "house_size", "project", "contractor", "deconMaterialName", "SimpleEOLname", "disposition", "quantity", "quantity_units", "material_impacts","site_transport_impact", "EOL_transport_impact", "impact_units"))

# join all 3 data frames to create a master data set that contains all material and transport impacts from both scenarios and save it
final_complete_data <- bind_rows(decon, demo, worker_equipment_impacts)
final_complete_data$project <- as.character(final_complete_data$project)
final_complete_data$contractor <- as.character(final_complete_data$contractor)

write.csv(final_complete_data, "intermediary/final_complete_data_all_scenario_impacts.csv")

# in some instances, it may be desirable to see material impacts of either just the Softwood lumber products, or everything else minus the Softwood lumber products
softwood_lumber <- final_complete_data %>%
   filter(deconMaterialName == "softwood lumber")

nonSoftwood_lumber <- anti_join(final_complete_data, softwood_lumber, by = c("scenario", "house_age", "house_size", "project", "contractor", "deconMaterialName", "SimpleEOLname", "disposition", "quantity", "quantity_units", "material_impacts", "site_transport_impact", "EOL_transport_impact", "impact_units"))

# remove comment to write/save these two tables separately as needed
# write.csv(softwood_lumber, "intermediary/softwood lumber data.csv")
# write.csv(nonSoftwood_lumber, "intermediary/nonSoftwood lumber data.csv")

```


The resulting `worker_equipment_impacts` data frame has `r length(worker_equipment_impacts$scenario)` rows and the same `r length(worker_equipment_impacts)` varible field columns as the `decon` and `demo` data tables. Now these three data tables (with exactly same column order and column names) can be combined into a single `final_complete_data` table and save as `final_complete_data_all_scenario_impacts.csv`. This final table contains `r length(final_complete_data$scenario)` rows of observations, and retains information specific to project, material weight quantity and material type category, disposition, impact category (designated by the `impact_units` variable), and impact values associated with the materials, site-transport (which includes worker and equipment per-project-transport impact values), and EOL-transport. This structure allows for the summary and break out of material based impacts, and all forms of transport impacts.


## Data Set Description

```{r data summary}
ave_house_size <- summarise(house_weight, avg_sq_ft = mean(house_size))
total_house_size <- summarise(house_weight, total_sq_ft = sum(house_size)) # total sq_Ft from all projects
```


**FIGURE X: Relation of disposal material weight to house size**

There was request from the project team to explore the weight ratio of salvaged materials to the dropbox disposal. This is depicted in FIGURE 3 where the dropbox total weight is on the x-axis, the salvaged material weight is on the y-axis, each project is a point along the x and y where the size of the point is relative to the size of the house, contractors are assigned to the point color, and each point is labeled with its respective percent salvaged material numerical value. In this situation, given the wide range and variation in the material weight values from each project, the median value is a more appropriate measure of central tendency to express an *average* value because it is not influenced by extreme values. For this reason, median values are included in the plot in order to provide context for the data points in terms of the percent salvaged materials and the total salvage and disposal weights. The median crosshairs have been added to visually assist in differentiating the better than average performers in the upper left quadrant (with higher than average salvage weight and lower than average dropbox weight), from the lower than average performers in the lower right quadrant (lower than average salvage weight and higher than average dropbox weight)

```{r ratio of salvage to dropbox weight}
ggplot(house_weight, aes(total_dropbox_quantity, total_salvage_quantity, color = contractor, size = house_size)) +
   geom_jitter() +
   #scale_size("house_size (sqft)") +
   geom_text_repel(aes(label = as.integer(percent_salvaged)), seed = 19, show.legend = F) +
   #geom_rug(aes(size = 0.25), sides = "bl", position = "jitter", show.legend = F) +
   geom_vline(xintercept = median(house_weight$total_dropbox_quantity), color = "grey30" ) +
   geom_hline(yintercept = median(house_weight$total_salvage_quantity), color = "grey30" ) +
   scale_radius("house size (sqft)", range = c(2,6)) +
   scale_color_calc() +
   labs(x = "kg dropbox materials", y = "kg salvaged materials") #+
   #ggtitle("Ratio of salvaged materials to dropbox materials by house size")

#ggsave(filename = "ratio_salvage_dropbox.png", device = "png", path = graph_path, units = "in")
# no need to save it here since the next document produces and saves the same graphic
```

**FIGURE _ : Salvaged weight to dropbox weight ratio**  



# Methods
The work process for the project is illustrated by the following figure. Data frames developed from the cleaned raw data summary (`01_decon_material_weight_conversion.Rmd`) and the data preparation notebook (`02_decon_data_prep2018.Rmd`) are used with the impact factors from LCA data and transport data collected from Oregon DEQ and Oregon Metro (`03_Transport.Rmd`) in the development of this data analysis notebook. This notebook provides a detailed account of the analysis and is used to produce custom R Markdown reports for different stakeholders based on the stakeholder information needs. 
  

![Project Workflow](data/project_workflow.png)

## Data Variable Distribution tests

Where we are interested in statistics of the variables in the data, the distributions of the data must be tested for assumptions before we can apply linear modeling methods. The Shapiro-Wilk's method for normality test is used on the `house_size`, `house_age`, and salvage and disposal weight variables.

```{r Shapiro-Wilk test}
shapiro.test(house_weight$house_age)
shapiro.test(house_weight$house_size)
shapiro.test(house_weight$total_salvage_quantity)
shapiro.test(house_weight$total_dropbox_quantity)
```

```{r quantile plots}
ggqqplot(house_weight$house_age) + ggtitle("house_age qqplot for normality")
ggqqplot(house_weight$house_size) + ggtitle("house_size qqplot for normality")
ggqqplot(house_weight$total_salvage_quantity) + ggtitle("total_salvage_quantity qqplot for normality")
ggqqplot(house_weight$total_dropbox_quantity) + ggtitle("total_dropbox_quantity qqplot for normality")
```





# Results


## Data Variable Relations and Linear Correlations

```{r lin reg scatter}
# add lm of dropbox total kg ~ house size w/ r^2
# house_weight$contractor <- as_factor(house_weight$contractor)
house_in_dropbox <- lm(total_dropbox_quantity ~ house_size, data = house_weight)
summary(house_in_dropbox)
#plot(house_in_dropbox, 1:6)

house_salvaged <- lm(total_salvage_quantity ~ house_size + house_age, data = house_weight)
summary(house_salvaged)
#plot(house_salvaged, 1:6)
```


```{r correlations}
num_house_weight <- house_weight %>%
   select("total_salvage_quantity", "total_dropbox_quantity", "house_age", "house_size")
names(num_house_weight) <- c("salvage kg", "dropbox kg", "house age", "house size")

house_weight_corr <- cor(num_house_weight, method = "pearson")
print(house_weight_corr)
cor.test(house_weight$total_salvage_quantity, house_weight$house_size, method = "pearson")
cor.test(house_weight$total_dropbox_quantity, house_weight$house_size, method = "pearson")

corrplot::corrplot(house_weight_corr, method = "color", type = "lower", addCoef.col = "black", tl.col = "black", tl.srt = 0, sig.level = 0.05, insig = "blank", diag = FALSE, mar = c(0,0,1,0))

```


House size is not a valid predictor of the total dropbox quantity since only a fraction of a percent of the variance in the house sizes explains variance in the dropbox quantity and the correlation is not statistically significant. A closer look at a correlation matrix between the salvage and dropbox quantities and the house characteristics reveals a modest positive correlation of 0.41 between the square foot size of the house and the salvage quantity, which when tested, turns out to be a statistically significant correlation at the 0.01 alpha level (a 0.05 p-value is conventionally accepted as an indicator of statistical significance). 



For the following net impacts calculations, each impact is summarised by impact category and activity. The sets of `net` impacts are obtained by subtracting the `decon` activity from the `demo` activity: $$demo - decon = net$$.

The `net` impacts are then expressed as a `net_per_sq_ft` value,  which is `net` divided by the sum total of square feet of the projects (`r sum(house_weight$house_size)`) and a `net_ave_project` value, which is `net` divided by the total number of projects (`r length(house_weight$project)`): 
$$\frac{demo - decon} {total ~~ sq ~~ ft} = \frac{net}{sq~ft}$$ 
$$\frac{demo - decon} {total ~~ number ~~ of ~~ projects} = average ~~ project ~~ net ~~ benefit$$

These calculations are completed for each impact category and saved as separate .csv files.

```{r worker equipment summary}
# for the final data summaries, Jordan will want to see this broken out from the site transport leg
# # create a table showing demo and decon worker and equipment total impacts and later add per sq_ft and per project

worker_equipment_summary <- worker_equipment_impacts %>%
   group_by(scenario, deconMaterialName, impact_units) %>%
   summarise(impact = sum(site_transport_impact)) %>%
   rename(impact_origin = deconMaterialName)
```


Reshape the `final_complete_data` table for use with graphics and table construction & final results development.

```{r gathered impacts}
gathered_impacts <- final_complete_data %>%
   gather(`material_impacts`, `site_transport_impact`, `EOL_transport_impact`, key = "impact_origin", value = "impact")
```

Calculate net impacts by impact category and save tables.

```{r net impact summary table}
net_impacts <- filter(gathered_impacts, SimpleEOLname != "per_project_transport") %>%
   group_by(scenario, impact_units, impact_origin) %>%
   summarise(impact = sum(impact)) %>%
   full_join(worker_equipment_summary, by = c("scenario", "impact_units", "impact_origin", "impact"))

# add the per house column now and export data table for the report draft
# per comments from Jordan, the next several code chunks showing net impacts will not be used,
# instead, make new tables showing impacts by scenario per average home with net value
final_impact_summary <- net_impacts %>%
   mutate(per_home_impact = impact / length(house_weight$project))

write.csv(final_impact_summary, "intermediary/final_impact_summary.csv")
```


Carbon net impact summary:
```{r carbon net impact summary}
carbon_net_impacts <- filter(net_impacts, impact_units == "kg CO2e") %>%
   ungroup() %>%
   spread(key = scenario, value = impact)

carbon_net_impacts <- carbon_net_impacts %>%
   add_row(impact_units = "kg CO2e",
           impact_origin = "total",
           decon = sum(carbon_net_impacts$decon),
           demo = sum(carbon_net_impacts$demo)) %>%
   mutate(net = demo - decon) %>%
   mutate(net_per_sqft = net/as.numeric(total_house_size)) %>%
   mutate(net_per_project = net/length(house_weight$project)) %>%
   add_column(order = c(3,5,1,2,4,6)) %>%
   arrange(order) %>%
   select("impact_origin", "decon",  "demo", "net", "net_per_sqft", "net_per_project", "impact_units") 
   
write.csv(carbon_net_impacts, "output/carbon_impact_breakdown.csv")

kable(carbon_net_impacts, caption = "Carbon Impact Table")
```


Energy net impact summary:
```{r energy net impact summary}
energy_net_impacts <- filter(net_impacts, impact_units == "MJ") %>%
   ungroup() %>%
   spread(key = scenario, value = impact)

energy_net_impacts <- energy_net_impacts %>%
   add_row(impact_units = "MJ",
           impact_origin = "total",
           decon = sum(energy_net_impacts$decon),
           demo = sum(energy_net_impacts$demo)) %>%
   mutate(net = demo - decon) %>%
   mutate(net_per_sqft = net/as.numeric(total_house_size)) %>%
   mutate(net_per_project = net/length(house_weight$project)) %>%
   add_column(order = c(3,5,1,2,4,6)) %>%
   arrange(order) %>%
   select("impact_origin", "decon",  "demo", "net", "net_per_sqft", "net_per_project", "impact_units") 
   
write.csv(energy_net_impacts, "output/energy_impact_breakdown.csv")

kable(energy_net_impacts, caption = "Energy Impact Table")

```


## Salvaged Materials Overview





```{r house weight, include = F}
ggplot(gathered_impacts, aes(SimpleEOLname, impact, fill = contractor) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_hline(yintercept = 0,  color = "grey0") +
   scale_fill_manual(values = DEQ_pal) +
   coord_flip() +
   facet_grid(scenario ~ impact_units) +
   ggtitle("impacts by contractor and basic material type")

ggplot(gathered_impacts, aes(project, impact, fill = SimpleEOLname) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_hline(yintercept = 0,  color = "grey0") +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   theme(axis.text.x = element_blank()) +
   facet_grid(scenario ~ impact_units) +
   ggtitle("impacts by project and basic material type")

ggplot(gathered_impacts, aes(project, impact, fill = SimpleEOLname) ) +
   geom_bar(stat = "identity", position = "fill") +
   geom_hline(yintercept = 0,  color = "grey0") +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   theme(axis.text.x = element_blank()) +
   facet_grid(scenario ~ impact_units) +
   ggtitle("impacts by project and basic material type")
```



## Transport Factors and Equipment Use

```{r impacts by scenario and material, include = F}
ggplot(gathered_impacts, aes(x = scenario, y = impact, fill = impact_origin) ) +
   geom_bar(stat = "identity", position = "dodge") +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   labs(x = "scenario", y = "impact total") +
   facet_grid(impact_units~.) +
   ggtitle("Impacts of materials and transport")

ggplot(final_complete_data, aes(x = scenario, y = material_impacts, fill = SimpleEOLname) ) +
   geom_bar(stat = "identity", position = "dodge") +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   labs(x = "scenario", y = "impact total") +
   facet_grid(impact_units~.) +
   ggtitle("Impacts by simple material category")
```



## Contribution Analysis
In order to show material flows and their impact contributions, the gathered data set needs to be summarised.  
```{r grouped & gathered data, include = F}

ggimpact <- gathered_impacts %>%
   group_by(scenario, impact_units, impact_origin, deconMaterialName, SimpleEOLname) %>%
   summarise(impact_sum = sum(impact)) %>%
   filter(impact_sum != 0 )

#write.csv(ggimpact, "intermediary/Impact summary by origin and material type.csv")

```

Then for each material, material type, impact category, and origin of impacts from materials and various transport factors, the sum of the impact values can be displayed graphically to show the specific levels of impact contribution to the total impact value.
```{r impact contributions, fig.height = 12, fig.width = 14}
ggplot(filter(ggimpact, impact_units == "MJ"), aes(x = deconMaterialName, y = impact_sum, fill = SimpleEOLname) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_text_repel(aes(label = as.integer(impact_sum) )) +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   coord_flip() +
   ggtitle("Energy impacts (MJ) by origin material type") +
   facet_grid(impact_origin ~ scenario)

ggplot(filter(ggimpact, impact_units == "kg CO2e"), aes(x = deconMaterialName, y = impact_sum, fill = SimpleEOLname) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_text_repel(aes(label = as.integer(impact_sum) )) +
   scale_fill_manual(values = DEQ_pal[5:9]) +
   coord_flip() +
   ggtitle("Carbon impacts (kg CO2e) by origin material type") +
   facet_grid(impact_origin ~ scenario)


# use the next two sets of ggplot calls to experiment with creating a graphical chart of all impact origins & material types
ggplot(filter(ggimpact, impact_units == "MJ"), aes(x = deconMaterialName, y = impact_sum, fill = impact_origin) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_text_repel(aes(label = as.integer(impact_sum) )) +
   scale_fill_manual(values = DEQ_pal[12:14]) +
   coord_flip() +
   ggtitle("Energy impacts (MJ) by origin material type") +
   facet_grid(SimpleEOLname ~ scenario )

ggplot(filter(ggimpact, impact_units == "kg CO2e"), aes(x = deconMaterialName, y = impact_sum, fill = impact_origin) ) +
   geom_bar(stat = "identity", position = "stack") +
   geom_text_repel(aes(label = as.integer(impact_sum) )) +
   scale_fill_manual(values = DEQ_pal[12:14]) +
   coord_flip() +
   ggtitle("Carbon impacts (kg CO2e) by origin material type") +
   facet_grid(SimpleEOLname ~ scenario)

```

